Multicomponent reaction-diffusion processes on complex networks 
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We study the reaction-diffusion process A + B — » on uncorrelated scale-free networks analytically. By a 
mean-field ansatz we derive analytical expressions for the particle pair-correlations and the particle density. 
Expressing the time evolution of the particle density in terms of the instantaneous particle pair-correlations, we 
determine analytically the 'jamming' effect which arises in the case of multicomponent, pair-wise reactions. 
Comparing the relevant terms within the differential equation for the particle density, we find that the 'jamming' 
effect diminishes in the long-time, low-density limit. This even holds true for the hubs of the network, despite 
that the hubs dynamically attract the particles. 

PACS numbers: 89.75.Hc, 82.20.-w, 05.40.-a 
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I. INTRODUCTION 

The merging of graph theory and methods from classical 
statistical physics has led to the modern theory of complex 
networks |ll|2j]. This development went along with the dis- 
covery of complex networks in a vast variety of sciences. Al- 
most all observed real networks share a so called scale-free 
degree distribution. That is, the number of links (the degree 
k) of a randomly chosen vertex of the network is distributed 
according to 



P(k) « k-\ 



(1) 



For exponents 2 < y < 3 this distribution has a diverging sec- 
ond moment (k 2 ) for infinite size networks (N — > °°). A very 
prominent example is the WWW |.3j with a value y of ap- 
proximately 2.1. As a consequence, the network is very ro- 
bust against the removal of random nodes, but extremely weak 
upon removal of the highly connected nodes |4j] . 

The diverging (k 2 ) leads to enormous fluctuations in the de- 
gree of a vertex and causes a large heterogeneity in the vertex 
connectivity. This crucial feature of scale-free networks has 
strong consequences on dynamical processes taking place on 
them. The absence of an epidemic threshold in the presence 
of a diverging (k 2 ) is one of the most notable examples JH^]. 

A similar class of dynamical processes, the reaction- 
diffusion proces ses B|, has been studied numerically 
on networks by iGallos and Argvrakisl |0, Q3] and shortly 
after also analytically within mean-field (MF) theory by 
ICatanzaro et alJ full (see fl2il for a recent series-expansion 
approach to networks). This type of dynamics is capable of 
modeling epidemic spreading, chemica l reactio ns, and many 
more. The analytical work by ICatanzaro et all calculated the 
density decay of the one component A + A — > reaction for 
homogeneous and heterogeneous complex networks. The re- 
sults for homogeneous networks agree with the classical MF 
behavior of a density decay linear in time, 1 /p^ (f ) « t . In the 
case of heterogeneous (i.e., scale-free) networks, one finds a 
more general power law for the density decay 
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with the y-dependent characteristic exponents given by 

a(y) = ( 1/(y - 2) 2<y<3 (3) 
y " [l y>3 

and 



P(Y) = 



y=3 
otherwise 



(4) 



pM(f) 



oc t 



a(y) 



(lnf) 



P'.y) 



(2) 



The reason for this extremely fast density decay originates in 
the existence of a few vertices with a very large degree (so- 
called hubs) in a scale-free network. As the analysis shows, 
the density is constantly 1 /2 for vertices with a degree k > 

* c = (t)/2pW(0. 

If one wants to extend these results to more complicated 
diffusion-annihilation dynamics, for instance the A+B — > 
process, which up to now has been studied only numerically 
191 1 1 0fl . one realizes that this diffusion-annihilation dynamics 
introduces forbidden diffusion steps, since a vertex can hold 
at most one particle at a time. This constraint prohibits an A 
(B) particle to diffuse to an adjacent vertex which is already 
occupied by another A (B) particle. As we are going to discuss 
below, this leads to a 'jamming' effect, a new feature of the 
dynamical process. Nevertheless, the obtained solution has a 
structure similar to that of the A + A — > process. The paper 
begins in Section|ll]with a MF calculation of the particle pair- 
correlations in dependence of the particle density. It follows 
a MF ansatz for the density decay of the A + B — > process 
from which we quantify the 'jamming' effect. The results of 
the analytical work are compared with numerical simulations 
in Section^] At last, we conclude with Section HVl 



II. THE A+B — » REACTION ON COMPLEX NETWORKS 

We assume the complex network of nodes to be fully 
defined by its N xN adjacency matrix an. To discuss a phys- 
ically meaningful complex network in the sense of diffusion, 
we take the network to be undirected and free of self- and 
multiple-connections. Therefore, aij is a traceless, symmetric 
binary matrix with the elements a, ; - being or 1, which sym- 
bols a (dis)connection between site i and j. The state of vertex 
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i at time t is described by two dichotomous variables nf\t) 

and nf\t). Their values may be 1 or only, indicating the 

presence or absence of a particle A («•(?)) and B («.(?)). 
The system state is thus defined by 



nT\t)){nf{t)), we obtain 



n(f) =n (a) (f)+n( b) (?) 
n^Xt) = {nf>{t),4\t),...,4\t)} 
nV{t) = {n^Xt),nf{t),..J*Xt)}. 



(5) 



Note, that these variables have to fulfill the constraint 
n- (f)n,- (f) — at any time. In the course of the calcula- 
tion, we will take the average over multiple realizations of the 

same system, turning the discrete nfXt) and nf'Xt) variables 

into densities p. (f ) and p| b ^ (f ). Furthermore, we will assume 
throughout the analysis the statistical equivalence of vertices 
of the same degree k. Therefore, denoting by V (k) the set of 
all vertices with the same degree k, we assume that 



pW(0 = p? } (0 Viev(t) 
p, (b) (f)^p[ b) (0 Viey(t) 



(6) 



is valid. Following standard MF treatment, we hence ne- 
glect all fluctuations which might exist within a set of vertices 
V (k). The total density p (a) (f) (p^ b ^(f)) is given by the set 

of partial particle densities {pj^(f)} ({p* (*)}) through the 
relation 



) w=£pf ) wn*)- 

k 



(7) 



A dynamics starts by random assignment of maximal one par- 
ticle per vertex. The particles diffuse by random jumps at a 
rate X to adjacent neighbors through the network. If two dif- 
ferent particles meet at a vertex, they instantly annihilate and 
the vertex becomes empty. Before we proceed to the time 
evolution of the system, we first derive an expression for the 
particle pair-correlations in dependence of the partial particle 
densities. 



A. Particle pair-correlations 

We quantify the particle pair-correlations for given partial 
particle densities by counting the number of contacts between 
particles on adjacent vertices. To count the AB contacts of a 
vertex i, we assume vertex i to carry an A particle and count 
the number of adjacent vertices which are occupied by a B 
particle. Setting this number in relation to all connections of 
the vertex i yields the pair-correlation coefficient 



<^\t) = U*Xt)W; 
1 j 



(b) 



(8) 



Averaging now over a whole ensemble of equal systems and 
making use of the usual MF assumption (nfXt) n.j(t)) ~ 



Qt\t) = T.P?Xt)I,aijP?Xt) 



(b), 



(9) 



By using the statistical equivalence of all Nk vertices i with the 
same degree k, we can sum over all these vertices such that 



ejf b) (0 = TP? ) Epjf ) (0^ I E 



(10) 



iev(k) jev(k') 



In this step, we split the sum with index j into two sums over 
k' and one over V(k'). The double sum over aij is related 
to the conditional probability P(k'\k) that a vertex of given 
degree k has a neighbor which has degree k' . This equation 
has been derived previously |6] to be 
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£ £ a u =P(k'\k). 



(11) 



iev(k) jev(k') 



Using this equation and assuming an uncorrelated net- 
work, which simplifies the conditional probability P(k'\k) to 
k'P(k') j (k), we obtain the expression 



e [ ab) (r) = p[ a) (r)0( b )( f ), 
where we define 

e (b) W = ^E*'P? ) WW 



(12) 



(13) 



One should note that by introducing the mean (k), the values 
of the exponent y are limited to y > 2. Otherwise the mean 
(k) is not defined in the limit of infinite size networks (N — > 
°°). The overall particle pair-correlation coefficient 2^ ab '(f) 
can easily be computed by multiplying Eq. d!2i with P(k) and 
summing once more over all k, 



Q(" b >(f)=p00(f)e0>)(f). 



(14) 



Analogously, we have gW(f) = p( a )(?)©( a )(?), 2< bb )(f) = 
p (b) (j ) ©(b) and qM (f) = p (b) (f ) @ (a) 



B. Density decay 

For further computations we will assume for simplicity that 
the initial densities of p( a ' and p( b ' are equal, so that there 
is a symmetry between A and B particles. We will calculate 
only an expression for n' a '(f), and one may obtain the corre- 
sponding «( b '(f) equations by interchanging indices A and B. 
Modeling the diffusion as a Poisson process 1 13], the set of 

fa) 

{n\ (f)} changes within an infinitesimal time interval df as 
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nr(t + dt) = n?>mt i \d t )+ l-( n ?>(t)+nr(t)) 



.00/ 



:(a), 



(15) 



Here T| ^ and f^f 1 are dichotomous random variables, taking values of or 1 with certain probabilities p and l—p respectively, 



(a) ( A \ J P = Xdt 

1\) >(dt) = 



a ijn f{t) ( l (a) \ 



1 l-p 



1 p = Xdf£ 

,0 l-p 



(a) 

fly/l) '(/) 



(16) 



(17) 



The following two cases need to be distinguished: (i) If site i is occupied by an A particle at instant t, r\) (df ) is responsible for 

the next time step: The site may become empty (r|j a ' = 0) with a probability proportional to the product of the jumping rate X 
and the time interval df if a B particle in the neighborhood jumps onto site i or if the A particle at / jumps away to a neighborhood 

(a) 

site where no A particle is already located. Otherwise no change happens, (ii) If the site i is empty at instant f, then q. (df) will 
determine the time evolution: The vertex may become occupied by an A particle only if one in the neighborhood jumps onto 

vertex i. Note that the two random variables T|} and are hence not independent from each other, but we will treat them as 
independent (cf. II 110 . 

fa) 

Equation 1151 yields an average time evolution for n) (t) 



E 

j 



(a), ^ijnf\t) (a) 1 (a) 



kj k 



nt\t)ja ijn f (t)-(l- «, laJ (f) + », (b) (f) 



cLijnJ (f) 



(18) 



where we have set without loss of generality the jumping rate X = 1 . Averaging over a whole set of equal initial configurations 
and applying once more Eq. il It and the statistical equivalence of vertices with the same degree, Eq. ©, we obtain after some 
formal rearrangements 



dp 



df 



~pf-E i 



(a) (b) , j (a) 1 (a) (a) . (a) (a) . (b) (a) 
Pi 'Pk' ~ k Pk lPk> ~Pk> +Pk Pk' +Pk Pk' 



kP{k'\k)\. 



(19) 



Here we have suppressed the explicit time dependence for the sake of simplicity. Assuming the network to be uncorrelated (i.e. 
that P(k'\k) = k 1 P(k') / (k)) allows us to perform the sum over k 1 , yielding finally the expression 



df 



k r 



Pi 



Wp(b)_p(a) + p WpW +p 0»)p(a)l +p W (.) 



(20) 



for the partial particle densities. Multiplying Eq. J20i with 
P(k) and summing over all k values results in the differential 
equation for the overall density, 



dp 



(a) 



df 



= _p(b) (a)_ p (a) (b) 



(21) 



Q(ab) _ q 



(ba) 



tides. To proceed further, we need expressions for p[ and 



Pk 
Pt 



(bi 



Since the initial densities are equal, we have forcibly 
= g(bb) because of symmetry. This implies the equality 
= p[ b ^ = p' k , allowing further simplifications and trans- 



From Eq. i2l\ it is apparent that the density decay is di- 
rectly proportional to the pair-correlations among unlike par- 
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forming Eq. i20\ into 

# = -Pl + |y[i-3pyp' + P^ 

= -p' k +^[l-3p> k ]p' + Ql 



(22) 



This differential equation is very similar to the one previously 
found for the A + A — > process lUlll . with an additional term 



Q'l 



Q)T' — an( i a coefficient of 3 instead of 2 in 

front of p[. The additional term measures the number of con- 
tacts among particles of the same type, which slow down the 
reaction as these jumps are prohibited, causing a 'jamming' 
effect. 

To test whether this new term Q'£ = p' k & alters the behavior 
of the dynamics, we compare it to the other term 3&p[p'/ (k), 
which is as well quadratic in the density, and find 



©' _ i L^'pW) _ i W> 



3p>k/(k) 3k ^p' k ,P{k') 3k (p^) 



(23) 



Due to the fact that the particles are dynamically attracted 
by the hubs of the network, which has been shown analyti- 
cally in Ref. llHl for the A + A — > process and numerically 
for the A + B — > process in Ref. |21 as well, we propose 
the following approximation to proceed: We know from the 
A + A — ► process that hubs (k > k c ) drive the dynamics, and 
that the density on those hubs is almost constant, p c = 1/2, 
whereas the density on vertices which are not hubs is substan- 
tially lower. If we approximate the densities of all vertices 
with a degree k < k c to be zero, the terms (k'p' k ,) and (p^,) 
become in the thermodynamic limit 

(kp' k ) « £ kP(k) « p' c Hk^dk = P' C $^ (24) 

k=k c Jk = ' 2 



l-Y 



<P*>«PcLn*)«Pc/ ™ = Pc^T- (25) 

k=kc Jk ° "~ 1 

In this step, we apply the continuous k approximation, which 
allows us to replace the sum by an integral. By inserting ex- 
pression J24t and J25l > into Eq. d23i . we find 



0' 



3p'k/(k) ~ 3y-2 k 



(26) 



Therefore, we can neglect Q'l for a vertex with k ^> k c . For 
k <C k c , we can neglect Q k (and 3kp' k p' / (k)) as this term is 
quadratic in the density, being very small for nodes with k <C 
k c in the long-time, low-density limit. The intermediate range 
of k w k c is difficult to assess analytically and needs to get 
quantified in the next Section by numerical simulations, which 
show that expression (I23> is substantially smaller than 1 in the 
low density limit even for the range where km k c . Therefore, 
'jamming' is only of relevance for vertices with a low degree 
and high densities. 

The calculation leading to Eq. ( I26> is carried out in the limit 
of infinite size network. In all real networks, one inevitably 
has a maximum degree £ max inducing finite-size effects. This 



maximum degree k max limits the upper bound of the integrals 
in Eqs. ( 1241 and fl25t . Evaluating these integrals with such an 
upper bound, one obtains 



0' 



„ iy-i*c 

3p'k/(k) ~ 3y-2 k 



f{k c /kmaxh 



with 



/(*) = 



i- x y- 1 ' 

The scaling function f(x) has the limiting values 



fix) = {1-2 
,Y-1 



(27) 



(28) 



(29) 



Since f(x) is a monotonically decreasing function for y > 2, 
the finite-size effect on the result in Eq. J26l > is to slightly de- 
crease the importance of the 'jamming' term evenly for all 
degrees k. 

Concluding that the 'jamming' effect is not relevant in the 
long-time, low-density limit, we neglect the Q k = p' k & term 
in Eq. d22i and obtain an equation which relates p' k to its 
derivative. However, we aim at a relation for p' k itself. As done 
for the A + A — > process 11 ill , we proceed with the quasi- 
static approximation, setting dp' k /dt w 0. This assumes that 
the diffusion process is at any time much faster than the anni- 
hilation reaction (see also 1 14]). This approximation should be 
valid in the case of low densities, when the complex network 
is sparsely populated and the number of diffusion events in 
a time interval df is much larger than the annihilation events, 
implying that the particles are always in an equilibrium state 
with respect to the degree distribution. Doing so, we get an 
approximate expression for pL 



Pk 



P'k/(k) 
\+3p'k/{k)' 



(30) 



This expression for p' k has the same structure as found for the 
A + A — > process except that the coefficient of p' in the de- 
nominator is 3 instead of 2. As the structure of the differential 
equations is the same as for the A + A — > process, we obtain 
the same scaling-behavior of Eqs. Q,® for each component. 
The new critical k c for which a vertex is sensed as a hub by the 
dynamics is k c = (k) /3p' = 2(k) /3p. Therefore, for the hubs 
in the system with k > k c , Eq. J30i is close to 1/3, which is 
completely consistent with MF, as this means that hubs are oc- 
cupied by approximately 1/3 of the time by each component 
A and B, and are empty for the remaining 1/3 of the time. 



III. NUMERICAL SIMULATIONS 

To test the analytical results obtained, we performed inten- 
sive numerical simulations of the A +B — » process on scale- 
free networks. The uncorrected, scale-free networks are gen- 
erated with the uncorrected configuration model (UCM) al- 
gorithm and have a size of Af = 10 6 if not stated 
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otherwise . The exponents y simulated are in the range of 2 . 1 to 
3.5, while we only present a suitable subset in the figures. In 
short, one draws for a network of size N a random number for 
each vertex according to the degree distribution P(k) °< k~Y, 
with an upper cut-off k m! , x = N 1 / 2 to ensure that the generated 
network is uncorrelated lUsll . The drawn number corresponds 
to the target degree of each node and can be understood as 
half-edges to be joint with other half-edges to form a connec- 
tion. This is done in the central loop, in which one draws 
randomly two half-edges and joins them if this neither creates 
a self-connection nor a multiple-connection. Upon a success- 
ful join of two half-edges, these two half-edges are dropped 
from the set of eligible half-edges. In any case, one continues 
with the central loop by drawing again two half-edges, and so 
forth. After the distribution of all half-edges, we only keep 
the largest component of the generated network. We use a 
minimum degree of k m \ n = 2 in our simulations, so that the 
largest component usually coincides with the full network. 
We have verified that the networks are indeed uncorrelated 
by obtaining the degree-degree correlation coefficient 
which has an absolute value smaller than 10~ 3 in all cases. 
On these networks, the dynamics is simulated in the follow- 
ing way: Initially a fraction p = 2p( a ) = 2p' b ' of randomly 
chosen vertices is selected, which we choose as p = 0.1 or 
0.95. Then, the algorithm assigns randomly an equal amount 
of A and B particles to the set of chosen vertices. After this 
initial setup, the diffusion-annihilation dynamics starts. First, 
a vertex which carries a particle and a random adjacent neigh- 
bor of this vertex are randomly selected. Three cases need to 
be distinguished: (i) If the neighbor vertex is empty, the parti- 
cle moves to the new vertex, leaving the initial vertex empty, 
(ii) If the neighbor vertex is occupied by a particle of the other 
type, an annihilation reaction occurs and both vertices become 
empty. Accordingly, the number of particles is decreased for 



,00. 



1. 



each particle type by one, »W — > n ( a ) — \ t and n' b ) 
(iii) If the neighbor vertex is occupied by a particle of the same 
type, then no jump occurs. In any case, the time is updated by 
t — * t + 1 /(«W +nC0), w here and «( b ) correspond to the 
values before the diffusion step, and one continues by select- 
ing randomly another vertex carrying a particle, and so forth. 

In order to obtain the system's typical behavior, we average 
over 50 independent dynamics on each graph and over 100 
independent graphs, making up a total of 5000 dynamics per 
data-point. 



A. Validation of the approximations 

To validate the analytical calculations developed in the last 
section, we have to verify the two central approximations 
made which are based on the assumption of a small particle 
density on the network. Furthermore, it is crucial to get an 
estimate which densities can be considered small enough for 
the validity of the approximations. Our first approximation 
was to neglect in Eq. J20i the 'jamming' term Q' k ' in compari- 
son to the other term quadratic in the density 3kp' k p' / {k) . We 
have shown the validity of this approximation analytically for 
vertices with a degree k 3> k c and k <^k c (the latter for low 




FIG. 1 : Plot of k c and k c of the numerically simulated A + B — > pro- 
cess for an exponent y= 2.75 and an initial density po = 0.95. Since 
fc c is increasing much faster than k c , 'jamming' becomes quickly ir- 
relevant. 
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FIG. 2: Ratio of k c and k c of the numerically simulated A + B — » 
process on networks of different exponents y as indicated, with (a) an 
initial particle density po = 0.95 and (b) po = 0.1. 



densities). To check the intermediate range k rj k c , we per- 
form numerical simulations. If we set the ratio in Eq. d23l 
equal to 1, we obtain a critical degree k c , 



3p'/(k) 



(31) 



which separates vertices whose 'jamming' term is less impor- 
tant than the other quadratic density term in Eq. (I22> from 
those vertices for which the 'jamming' term is at least of equal 
importance. The time-evolution of the particle density on 
vertices with a degree k 3> k c is not affected by 'jamming', 
whereas vertices with a degree of the order of k c or lower 
are affected. On the other hand vertices with a small degree 
do not contribute to the overall particle density at later times, 
since the hubs dynamically attract the particles and carry the 
highest density p'. = 1/3. Considering only vertices as hubs 
which have a degree k> k c = (k) /3p' , we have as a condition 
for 'jamming' not being relevant 



k c <g:k c 



(32) 



If condition j32t is fulfilled, there are no vertices left in the 
network which do carry a sufficiently high density and whose 
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Rearranging Eq. d30t into 



(33) 



FIG. 3: Numerical validation of the quasi-static approximation ac- 
cording to Eq. 1331 for the A + B — > process, exemplified by an 
exponent y= 2.5 and an initial density po = 0.1. The fc-classes are 
logarithmically joined, where class 1 corresponds to vertices of de- 
gree k = 2, class 6 to 8 < k < 10, class 1 1 to 33 < k < 42, class 16 
to 134 <k< 178, and class 21 to 563 < k < 750. A data-collapse is 
observed for t > 50. 



'jamming' term is important for the time-evolution of their 
pk(t)- In Fig. [flwe exemplified this condition for an expo- 
nent y = 2.75 and an initial particle density po = 0.95. Note 
that the curves are only drawn until k c reaches the value of 
the maximum degree k max present in the network. Once the 
value of k c is much greater than k c the 'jamming' effect is of 
no more relevance for the time evolution of the process for all 
vertices in the network, including those with k~k c . It is cru- 
cial to note that k Q grows much faster than k c in the course of 
the process. Therefore, the 'jamming' is continuously dimin- 
ishing during the dynamics and only important for low degree 
vertices carrying a high density in the beginning of the pro- 
cess. An equivalent criterion to test weather 'jamming' is not 
relevant is to check if the ratio of k c /k c is substantially smaller 
than 1. In Fig. Eta) we illustrate this for an initial density 
po = 0.95. Again, the individual curves are only drawn until 
k c reaches £ max . They start with a maximum value of almost 1, 
indicating the presence of 'jamming' and drop quite quickly 
well below 1. In Fig.|2jb) we show the same simulations start- 
ing but with a much smaller initial density of po = 0.1. Most 
importantly, these curves already begin at values well below 
1 and therefore there is never 'jamming' present in the dy- 
namics. The interesting intermediate increase of k c /k c for the 
initial density po = 0.1 (Fig.|2jb)) comes from the fact the dy- 
namical hubs start with a density p = 0.05 which is smaller 
than their long-time density p ! c = 1 /3. Therefore, all vertices 
with p' k < 1/3 and a degree k > k c will have increasing parti- 
cle densities p' k which enter ©' in Eq. ( 13 It . Once the dynamics 
has reached its long-time behavior, there are no more p[-terms 
in ©' which increase in magnitude, since the dynamical hubs 
carry the highest density in the network. 

The second approximation made to obtain an expression for 
p^, the quasi-static assumption dp^/df w 0, yielded Eq. J301 . 



leads to an expression where the right hand side (and conse- 
quently the left hand side as well) is independent of k if the 
approximation is indeed valid. Plotting the left hand side of 
Eq. < I33I for a couple different degrees k should yield a data- 
collapse onto a single curve. Fig. [3] illustrates this in the case 
of an exponent y = 2.5. The curves join quite nicely at roughly 
t « 50. Similar time points are obtained for other exponents 
y. We can therefore expect that the quasi-static approximation 
holds after this time. 



B. Density decay and pair-correlations 

The verification of the scaling behavior for each compo- 
nent as predicted by Eq. l|3} and Eq. @ turns out to be a hard 
numerical problem, as finite-size effects occur quite early. A 
detailed discussion and derivation of finite-size effects for the 
A + A — > process can be found in Ref . full . The typical 
'scale-free' behavior of the dynamics corresponds to a density 
de cay as a power-la w of the time with an exponent larger than 
1. ICatanzaro et alJ full showed that the characteristic 'scale- 
free' behavior of the dynamics is driven by the 'dynamical' 
hubs of the system, where 'dynamical' hubs are those vertices 
which have a degree k > k c °= l/p( a ). As the density is de- 
creasing monotonically in time, the number of 'dynamical' 
hubs decreases as well. Thus, if in a finite network of size 
N there are no 'dynamical' hubs left, the density decay turns 
over to a decay linear in time. The vertex with the largest 
de gree k max therefore limits the duration until the cross-over 
from 'scale-free' to non-' scale-free' behavior happens. On 
the other hand, the largest degree k max being possible for un- 
corrected scale-free networks which contain neither self- nor 
multiple-connections scales with the square root of the system 
size, fc max °< N x l 2 , making very large system sizes necessary. 

One might expect that the scaling exponent derived in MF 
theory for each component's density decay, a(y) = l/(y — 2), 
is missing some effects which slow down the reaction for ex- 
ponents y close to 2. Otherwise, the diverging a(y) for y — > 2 
would result in a diverging reaction speed. Forcibly, we expect 
to recover the scaling law given by oc(y) for y — > 3. Choosing 
a value of y which is smaller than but close to 3 has the con- 
venient side effect that the density decay is relatively slow, 
such that the dynamics will show a 'scale-free' behavior for a 
longer period of time with the appropriate density decay ex- 
ponent larger than 1 . We verify this assumption by simulating 
the A + B — > process for various exponents y while keeping 
the system size constant at = 10 6 . In Fig. 0]we show the 
resulting density decays. The curves deviate from a linear in 
time decay only for very short durations, but the amount of 
time with which each process behaves 'scale-free' increases 
with increasing exponent y. To recover the scaling behavior of 
p( a ), we chose an exponent y = 2.75 which corresponds to a 
value a = 4/3. In Fig. [5] we present the results for the density 
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FIG. 4: Density decay of the numerically simulated A + B — » pro- 
cess for an initial density po = 0.1 and different exponents y. With 
increasing exponent y, the density decay behaves 'scale-free' typical 
for a longer period of time. The solid line has a slope of 1 and is 
shown as a guide to the eye. 
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FIG. 5: Density decay of the numerically simulated A + B — > 
process for networks of various sizes, exemplified by an exponent 
y= 2.75 and an initial density po = 0. 1. The plot illustrates the strong 
finite-size effects of the dynamics. For y = 2.75 the process is suffi- 
ciently slow to show the theoretic slope of 4/3 (drawn as a solid line) 
for large networks with size N = 10 7 for two decades. 



decay p( a ' for various network sizes of up to N = 10 7 . Never- 
theless, even for a system size of N — 10 7 , the density decay 
shows 'scale-free' behavior only for less than two decades. 

A far more sensitive alternative check of the analytical cal- 
culations is given by the direct comparison of the analyti- 
cally calculated particle pair-correlations, Eq. 1141 . and the 
measured ones from numerical simulations, as these particle 
pair-correlations directly drive the annihilation rate. There- 
fore, we evaluate the ratio between the analytical predictions 
of Eq. J 1 4i and the corresponding numerically obtained val- 
ues for the pair-correlations. It is crucial to note that Q^ ab '(0 
can be understood as a function of all partial particle densi- 
ties p^\t) and p[ b ^(f), so that it depends only implicitly via 
these densities on time t. Approximating these partial particle 
densities via the quasi-static assumption, one may write the 
particle pair-correlations as a function of the inverse particle 
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FIG. 6: Ratio between the analytically calculated correlations Q^l 

and the numerically obtained one Qj£L , plotted as a function of the 
inverse density p for an initial density po = 0. 1 and different expo- 
nents y. The range of the inverse density shown is limited by the va- 
lidity of the quasi-static approximation, which is valid at sufficiently 
low densities, p _1 > 5 X 10 . However, the inverse density is also 
limited from above to p _1 < 10 3 , as the pair-correlation is quadratic 
in the particle density. Such observables can only be meaningfully 
resolved as long as p > A' -1 / 2 , which yields the condition p > 10~ 3 
for the given system size of N = 10 6 . 

density p _1 by using Eq. QUI. 

o^^W)V w^mw) p(n (34) 

The inverse particle density p _1 can be regarded as an alterna- 
tive measure of time, since p _1 is a monotonically increasing 
function in the course of the dynamics. In Fig.|6j we present 
the resulting ratios for various exponents y. The expected 
value of 1 for the ratio, which would mean perfect agreement 
of analytical and numerical particle pair-correlations, is quite 
well achieved for exponents y close to 3. Smaller exponents 
y yield ratios somewhat larger than 1, which however depend 
only weakly on p _1 . The reason for the larger ratios is pre- 
sumably that these networks have larger fluctuations in the 
network connectivity structure. These topological fluctuations 
may lead to strong density fluctuations which are not captured 
by the current MF ansatz. For y > 3, the ratios become some- 
what larger than 1 as well, as it is seen for y = 3.5. This be- 
havior results from a segregation of the components, which 
is well known for lattice s |8] and has already been observed 
by Gal los and Argvrakia |9|] for networks with y > 3. That is, 
for networks with y > 3 the dynamics behaves as on lattices, 
which is as well reflected by the scaling relations from Eq. Q 
and Eq. 10} giving a linear in time density decay for networks 
with y > 3. Such segregation is not captured in our MF calcu- 
lation of the pair-correlations and is therefore out of range of 
its validity. 

IV. CONCLUSIONS 

This paper presents a detailed discussion of the A + B — > 
process on networks. We show analytically the existence of a 
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'jamming' effect for this two component reaction and quan- 
tify it analytically with the correlations among unlike parti- 
cles. The analysis can easily be generalized to multicompo- 
nent, pairwise processes and be used to derive a MF theory for 
a process where n particle types react pairwise with each other. 
Note, however, that the applicability is limited to cases where 
it is guaranteed that the particle densities are of the same or- 
der for all times. Otherwise, one component of the system 
with the largest density might cause a non-negligible 'jam- 
ming' and could significantly slow down the reaction. For the 



A +B — > process discussed here, we show that 'jamming' is 
only important for vertices with small degree k and high den- 
sities. In the long-time, low-density limit we derive analyti- 
cally that 'jamming' vanishes for all vertices in the network, 
including the hubs. This conclusion is supported by numerical 
simulations and allows us to reason that the particle densities 
of the A + B — > process on scale-free, uncorrected network 
show for each component the same scaling behavior as the one 
of the A + A — > process. 
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